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ABSTRACT 

This paper presents a hydrodynamical model describing the evolution of the 
gas reinserted by stars within a rotating young nuclear star cluster (NSC). We 
explicitly consider the impact of the stellar component to the flow by means of 
a uniform insertion of mass and energy within the stellar cluster. The model 
includes the gravity force of the stellar component and a central supermassive 
black hole (SMBH), and accounts for the heating from the central source of 
radiation and the radiative cooling of the thermalized gas. By using a set of 
parameters typical for NSCs and SMBHs in Seyfert galaxies our simulations show 
that a filamentary/clumpy structure is formed in the inner part of the cluster. 
This "torus" is Compton thick and covers a large fraction of the sky (as seen from 
the SMBH). In the outer parts of the cluster a powerful wind is produced, that 
inhibits the infall of matter from larger scales and thus the NSC-SMBH interplay 
occurs in isolation. 

Subject headings: AGN galaxies: nuclear starburst — black holes — hydrody- 
namics: accretion — methods: numerical 

1. Introduction 



The origin of the obscuring matter in active galactic nuclei (AGN) is one of the main 
challenges in modern astrophysics. The Unified scheme for AGNs requires of a dusty torus to 
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explain the two main types: non-obscured (type 1) and obscured (type 2) AGN (Antonucci 



1993 Urry & Padovani 1995). The model claims that the highly energetic central activity 



of AGNs is powered by mass accretion onto a supermassive black hole (SMBH) with mass 
MsMBR = 10^ — 10^^ Mq, and that the distinction between different types of AGNs is simply 
determined by the orientation of the torus. The origin, structure, and dynamics of the torus 
remains as one of the key unsolved problems in AGN physics, which presumably is also related 
to the SMBH feeding and feedback. Some models assume uniform gas distribution for the 
obscuring torus, whose thickness is supported by IR radiation pressure (e.g. Pier & Krolik 



1992 Krolik 2007 and references therein). Other models assume that the circumnuclear 



obscuring medium is clumpy, see for example Krolik & Begelman ( |1988 ) and Nenkova et al. 
(2002). However, the mechanism supporting the vertical thickness (in homogeneous models) 



and motions of the clumps (in clumpy models) is still under debate. Elitzur & Shlosman 



(2006) suggested that the torus is simply an outflow of dusty and optically thick clumps 



coming from the accretion disc . Konigl & Kartje (1994) proposed a model where a magneto 



centrifugal wind is responsible for the obscuration. Dorodnitsyn et al. (2011) proposed a 



model in which the obscuration is produced by a dusty wind driven by infrared radiation 
pressure from a dense torus. The torus as a wind model does not suffer from the vertical 
structure problem, but the origin of the wind is still unclear. Nayakshin et al. ( |2012 ) suggest 
that the obscuring matter comes from fragmentation of solid bodies (asteroids, comets and 



terrestrial- like planets) in the vicinity of the SMBH. Wada (2012 ), based on three dimensional 



simulations including radiative feedback from the AGN, describe the formation of a turbulent 
torus from the interaction of back-flows in a bipolar fountain, starting from a preexisting 



rotationally supported thin disk. Li et al. (2012) provided numerical simulations for the 



accretion flows with angular momentum which is sufficient to inhibit the accretion if the 
viscous processes are negligible, and to form a torus from the centrifugal gas. Most of the 
above models are dedicated mainly to explain the mechanism for vertical support of the 
torus and they are based on the assumption that a cold gaseous structure already exists. 
They also neglect the effects that the stellar feedback may provide to the inflow onto the 
SMBH. 



Wada & Norman (2002) and Schartmann et al. (2009) suggested that clumpy tori form 



with gas reinserted by stars within massive nuclear clusters during the late (post-supernovae) 
evolution. This model was motivated by the fact that stellar activity in the vicinity of 
the central SMBH has been found in a variety of Hubble type galaxies (see, for example, 
Fihppenko Olol[2003t [Gonzalez Delgado et all|2008| |Seth et al.|[2008| [Kormendy fc Bender 



2009, and references therein). Such nuclear starbursts are compact (< 50 pc) and have masses 
in the range 10^ - 10^ M© ( |Davies et all2007| [Watabe et aLl|2008t |Seth et al.|[20Tol ). Some of 



them show complicated star formation histories (Walcher et al. 2006) and present evidence 
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for global rotation, up to 45 km s ^, and stellar velocity dispersion of up to 120 km s ^. The 



Schartmann et al. (2009) model seems to be in agreement with Davies et al. (2007) and Wild 



et al. (2010) who claim that the accretion rates and thus the AGN luminosities rise rapidly 



at the late stage of the nuclear starburst evolution. However, this model does not explain the 
coincidence of luminous AGNs with young (ages less than 40 Myr) nuclear starbursts, as it is 
the case of NGC 1097 (see Figure 11 in Davies et al. 2007). Here we show that massive and 
compact young NSCs with a central SMBH can form filamentary/clumpy gaseous tori, and 
that the sizes of such tori are in the range inferred by near-infrared observations of Seyfert 



galaxies (Jaffe et al. 2004 Prieto et al. 2005) and with the sizes derived from their spectral 



energy distributions (SEDs, Alonso-Herrero et al. 2011, and references therein). These tori 



are consistent also with the compact sizes of the Compton thick medium estimated from 



X-ray observations of Seyfert 2 galaxies, type 2 AGNs (Risaliti et al. 1999). 



The paper is organized as follows. Section 2 deals with all the ingredients of the hydro- 
dynamical model. Details regarding the numerical scheme are given in section 3. Section 
4 describes the hydrodynamic solution which leads to the formation of the torus. We start 
the discussion with a model without the central radiation field. Then we include the cen- 
tral source of radiation and calculate column densities and possible obscuration fraction of 
the torus. We also compare analytic estimates of the centrifugal barrier with the numerical 
results. Section 5 discusses the impact of the NSC wind on the host galaxy and section 6 
presents our conclusions. 



2. The physical model 



We consider the accretion flow onto a SMBH embedded into a young stellar cluster 
which rotates like a solid body. The key point for our model is that the matter reinserted 
within young, massive, and compact star clusters could evolve in a catastrophic cooling 



regime which is very different from the Chevalier & Clegg (1985) adiabatic solution (see 



Sihch et al.|[2004| |Tenorio-Tagle et al.| [20071 [Wiinsch et aL][2008| hereafter S04, TT07 and 
W08, respectively). In the Chevalier and Clegg solution the thermalization of the matter 
ejected by massive stars inside the cluster leads to a high central pressure with an outward 
pressure gradient that steadily accelerates the gas from zero velocity at the center (i.e., the 
stagnation point is at the center) to the sound speed at the cluster edge. The reinserted 
matter exits the cluster as a free wind approaching its terminal velocity, which is twice the 
sound speed at the cluster border. If radiative cooling is considered, the gas in the central 
zones of massive clusters cools down and eventually becomes thermally unstable. This is 
because the average density of the gas increases linearly with the cluster mass while the 
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cooling rate inside the cluster volume grows as a square function of the stellar cluster mass. 
Consequently, the central pressure drops and the gas cannot be accelerated outwards. The 
stagnation point moves out of the cluster center and the solution becomes bimodal: within 
the stagnation volume the thermal instability leads to mass accumulation, while in the outer 
parts of the cluster a stationary cluster wind is established. The size of the stagnation zone 
becomes larger as one considers more massive clusters where strong radiative energy losses 
favor the frequent generation of cold parcels of gas (see TT07 for the semi-analytic procedure 
for calculating the stagnation radius). This leads to mass accumulation and eventually to 
star formation and thus to a positive star formation feedback (see |Tenorio-Tagle et al.|[2005| . 



In the presence of a central SMBH the stagnation point is always out from the center 



(Silich et al. 2008 Hueyotl-Zahuantitla et al. 2010, here after SOS and HZIO, respectively) 



and thus the flow in the vicinity of the central SMBH is always bimodal. Matter inserted 
within the stagnation volume forms an accretion flow whereas the mass inserted between the 
stagnation radius and the cluster edge drives the cluster wind. In such cases, the stagnation 
radius, Rst, is defined by the balance between the outward pressure gradient (strongly affected 
by radiative cooling) and the gravity force that makes the accretion flow very different from 
the classic Bondi solution. All matter reinserted within the stagnation zone remains bound 
to the cluster and thus defines the upper limit to the mass accretion rate onto the SMBH. 
In thermally unstable bimodal radiative cooling becomes more important, the 

stagnation point moves to a larger radius increasing substantially the mass accretion rate. 
Here by using typical values for the masses and sizes of NSCs, and masses of the SMBH 
in Seyfert galaxies, we explore the formation of the torus from the mass inserted within a 
rotating young NSCs with a central SMBH. 

Our physical model consists of a spherically-symmetric young NSC of radius -Rnsc and 
mass Mnsc with a homogeneous distribution of stars and with a central SMBH of mass 
-^SMBH- It accounts for the gravity pull from the stellar component and from the black hole. 
We consider a constant injection of mechanical energy (Lnsc) and mass (Mnsc) within the 
cluster volume via SNe II and stellar winds. Matter is inserted within the star cluster with 
a finite angular momentum given as the solid-body rotation around the polar axis: Vj-ot = 
ursm6(l), where u is the angular frequency, r and 6 are the radial and polar coordinates, 
respectively, and is the unit vector in the direction (p. We assume that the rotation velocity 
of the star cluster is small compared to the dispersion velocity of individual stars and thus 
we disregard the star cluster flattening due to its rotation. Radiative cooling is one of the 
main ingredients of the model and it is considered in all the computational domain. 

One of the main features of AGNs is their strong emission of ionizing radiation. A 
proper treatment of the emission from the central source requires to know the intensity 
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and the spectral energy distribution of the radiation field, and how it propagates through 
the ambient medium. Here we consider only the X-ray luminosity Lx, since such photons 
can penetrate deeper into high density regions shielded from UV photons. We consider a 
constant Eddington ratio Lx/L^dd = 0.08 during the calculations (LEdd ~ 1-3 x 10^^ erg 
s~^). This value is in the range 0.01-0.1 used by Wada (2012). In this way, the model 



includes Compton and X-ray heating from the central source and the radial component of 
the acceleration due to radiation pressure. These quantities were explicitly calculated using 
the ray-tracing method described in Appendix |A] 

The model implicitly accounts for shock heating by assuming full thermalization of 
stellar winds and SNe kinetic energy due to random interactions of the ejecta from massive 
stars that leads to gas temperatures of a few 10^ K. The reinserted gas at such temperature 
and relatively high density is not in thermal equilibrium. 

In all calculations we assume a maximum rotation velocity along the equator frot= 
50 km s~^ at the star cluster edge. The energy and mass deposition rates (Lnsc and 
-^Nsc) relate to the wind adiabatic terminal speed (Va,oo) by /^nsc = 0.5Mnsc^a,oo ~ 
3 X 1O'"'(Mnsc/1O^M0). The second expression arises when one scales the average results 



from Starburst99 (SB99, Leitherer et al. 1999). The wind adiabatic terminal speed is an in- 



put parameter in the model, we assume in all cases that it is constant and equal to 1000 km 
s~^. Note that this value is 2.5 times lower than the average adiabatic wind terminal speed 
(see 



Wiinsch et al.|[20TT ) for an instantaneous starbursts during the first 40 Myr. Therefore, 



Va,oo = 1000 km s ^ implicitly means additional mass loading to the fiow, which may be due 
to evaporation and destruction of preexisting high density molecular clouds and filaments. 



and/or evaporation of circumstellar discs forming low mass stars, see for example Stevens 



& Hartwell (2003); Melioli & de Gouveia Dal Pino (2004). We assume that mass loading is 



five times that inserted by massive stars. 



3. The numerical approach 



The numerical models presented here are based on the finite-difference Eulerian hydro- 



dynamic code ZEUS-3D version 3.5 (Stone & Norman 1992 Clarke 2010). All calculations 
were performed in spherical coordinates in 2D, with symmetry along the (p- direction. The 
set of hydrodynamic equations is 



dp 
di 



+ V ■ (pu) = g„ 



du 



+ (u- V)u = — VP 

P 



V$ + grad, 



'1^ 



(2) 
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de 

— + V ■ [u(e + P + p^)] =qe-Q + i/AGN, (3) 

and it is closed by the equation of state P = (7 — l)e, where e is the internal energy 
density and 7 = 5/3 the adiabatic index. The total internal energy is e = pv? /2 + e. The 
mass and energy deposition rates per unit volume are = 3MNsc/(4vrP^sc(l + ^71111)) and 
qe = 3LNsc/(47rP^gQ), respectively. The parameter r^mi represents mass loading. A small 
value of Va,oo means that the mass in winds of individual stars is loaded by additional 
mass from the parental cluster. The magnitude of the local acceleration due to gravity is 
I V$ 1= fifgrav = —GM{r)/r^, where M(r) = Msmbh + M^sci^ / RnscT is the mass enclosed 
within a sphere of radius r. The terms grad and Hagn in equations ^ and (|3|, respectively, 
represent the radial acceleration due to radiation pressure and the heating rate per unit 
volume due to the X-rays from the central source, the method used to calculate these terms 
is described in Appendix [A} The cooling rate per unit volume is Q = n'^A{T,Z), where 
n = p/p(T)mii is the gas number density, mn is the proton mass, and A(T, Z) is the cooling 
function, which depends on temperature T and metallicity Z. In all cases solar metallicity 
is assumed. To compute n we use an approximate treatment for the ionization degree: we 
take the values p(T > 10^) = 14/23 and p(T < 10*) = 14/11 as the mean mass per particle 
for ionized and neutral gas respectively. 

The model accounts for extremely fast cooling in all the computational domain and is 
considered in the calculation of the time step. The cooling rate is computed with an updated 
routine of W08 that uses, for temperatures above the ionization temperature of hydrogen. 



the Raymond & Cox coohng function tabulated by Plewa (1995), and for T< 10 K the 



Koyama fc Inutsuka| p)02[ ) cooling function: A(T) = Tki [10^exp(-1.184 x 10V(r+ 10^)) + 
1.4 X 10~^T^/^ exp(— 92/T)] erg cm^ s~^, where Fki = 2.0 x 10~^^ erg s~^ The minimum 
allowed temperature in the simulations is assumed to be T=100 K. Note that according to 



Joung & Mac Low (2006) the gas is thermally stable at temperatures where the slope of 



the cooling curve in the space log(A(T)) t>s T is > 1, what led Schartmann et al. (2009) to 
identify five stable zones in the Plewa (1995) cooling function. 

The flow was modeled following the prescription given in W08, which explicitly considers 
a continuous replenishment of mass and internal energy in all cells within the starburst 
volume at rates qm and qe, respectively. The inserted mass is subject to the gravity force of 
the SMBH as well as from the NSC, like in HZIO, however here the mass is injected with an 
angular momentum that corresponds to the solid body rotation of the cluster. In summary, 
the procedure applied to each cell within the cluster volume at every time step is: 

1. The radial velocity of the flow Vr is updated according to Vr = f r + (s'grav + gra.d)dt. 
where (^grav and (^rad are the local acceleration due to gravity and radiation pressure. 
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2. The density and total energy in a given cell are saved to poid and etot.ow- 

3. The mass is inserted so that pnew = PoW+^P, where 5p = {l+A^oiseOlmdi is the injected 
mass per unit volume. The mass 6p is inserted with rotation velocity frot, along the 
0— direction, assuming solid-body rotation for the star cluster: v^ot = ursinO, where 
u is the angular frequency, r is the distance from the center, and 6 the angle from the 
polar axis. 

4. The velocity is corrected so that the momentum is conserved: Vncw = VoidPoid/Pnew + 
0'i^rot'^p/Pnew, where the components of the velocity vector v are Vr, vg and v^. is the 
unit vector in the direction 0. 

5. The internal energy is corrected to conserve the total energy Cj^mid = etot.oid— PncwVn(,w/2- 

6. The new energy is inserted in a form of internal energy Cj^new = ej,mid + (l + ^noiseC)Q'e'^^- 

7. The AGN-heating i^AGN is included by increasing the internal energy in each grid cell 
by HAGNdt. 

In steps 3 and 6, C is a random number from the interval (-1,1) generated each time step 
it is used, and A^oise is the relative amplitude of the noise. The inclusion of the noise is 
necessary to break the spherical symmetry imposed by the initial conditions, an analysis of 
the effects of such noise is presented in W08 in the case of star cluster winds. Here we used 
their recommended value A^oisp = 0.1 in all our simulations. 



3.1. Initial and boundary conditions 

After an initial relaxation period, the solution reaches a steady state with a quasi- 
stationary wind blowing from the outer parts of the cluster and mass accumulation at an 
approximately constant rate in the inner region. To make the transition as short as possible, 
the models start from the stationary adiabatic wind solution of a star cluster with mass 
-^Nsc = IO^Mq, adiabatic wind terminal speed Va,oo = 1000 km s~^, and with radius in 
each case as given in Table [TJ The boundary conditions are set open at both r-boundaries 
and reflecting at both 6'-boundaries. We used the scaled grid option in r and a uniform one 
in 6. The computational domain and the number of zones in each direction were selected 
such that Ar ~ rA6, which conserves the shape of the zones and provides a higher resolution 
closer to the center. 
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Table 1:: Selected models and results from the simulations. 



No. 


^NSC 


MsMBR 




^rot 


Rad. 


Rst 


-Rx.num 




X24 


Xcold,22 




(Mo) 


(lO^Mo) 


(pc) 


(km s~^) 






(pc) 








1 


3.3 X 10*^ 


1 


10 


50 


N 


9.2 


~2 


1.95 


0.76 


0.92 


2 


3.3 X 10*^ 


1 


10 


50 


Y 


9.2 


~2 


1.95 


0.93 


0.86 


3 


3.3 X 10^ 


1 


10 


50 


Y 












4 


3.3 X 10*^ 


1 


40 


50 


Y 


31.3 


10.5 


9.7 


0.71 





5 


3.3 X 10^ 


10 


40 


50 


Y 













Note. — Summary of the models. The first cohimn marks the label of each model, where model 2 is our 
reference model; A^nsc and Msmbh are the masses of the NSC and the SMBH respectively; i?NSC is the 
radius of the star cluster; v^^t is the maximum rotation velocity at the star cluster surface; the labels Y and 
N, indicate whether a model includes the radiation from the central source or not. Rst is the 0— averaged 
value of the stagnation radius measured in the simulations at a time when the solution is steady; Rt is the 
radius of the centrifugal barrier, where the subscripts num and anl denote the results from the numerical 
simulations and the analytic prediction (Eq. [t]) respectively, numerically it corresponds to the time averaged 
radius of the point of maximum density; X24 indicates the time averaged of the fraction of the sky covered 
by column densities larger than 10^'' cm^^, and Xcoid,22 gives the fraction of the sky covered only by cold 
gas (T < 1500 K) with column density larger than 10^^ cm~^. 



4. Results 

We have calculated a set of models (see Table [T]) with NSC and SMBH parameters 



in the range given by the observations (see for example Figure 19 in Seth et al. (2010)). 
Models 1 and 2 have a star cluster radius -Rnsc= 10 pc and mass Mnsc = 3.3 x IO^Mq, 
respectively, both with a lO^M© SMBH. In model 1, the radiation from the central source 
is not considered. Hereafter we refer to model 2 as the reference model. In model 3 a star 
cluster is 10 times less massive than in the first two cases. Models 4 and 5 consider more 
extended clusters with i?Nsc= 40 pc, Mnsc = 3.3 x IO^Mq and SMBHs of Msmbh = IO^Mq 
and Msmbh = IO^Mq, respectively. The computational domain extends radially from 0.5 
pc to 15 pc in the case of models 1-3, and from 1 pc to 60 pc for models 4 and 5. The 
axial extent goes from 7r/2 — 1.3 to it/2 + 1.3 radians in all calculations. The complete set 
of relevant parameters for the simulated models are given in Table [1} 
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4.1. The hydrodynamic solution 

4-1.1. Case without a central source of radiation (model 1) 

We start with a description of model 1, in which the radiation from the central source 
is not considered. According to TT07 and W08 such NSC evolves in a catastrophic cooling 
regime. Here it is shown that the inclusion of gravity of the NSC+SMBH and the angular 
momentum of the inserted matter lead to the formation of a filamentary/clumpy torus. 




-15 -10 -5 5 10 15 -15 -10 -5 5 10 15 

z [pcj z [pc] 



Fig. 1. — : Filamentary torus resulting from model 1, snapshots at 1 Myr. The torus is 
composed by a collection of cold filaments and a dense core located at the centrifugal barrier, 
at 2 pc. In the left column, from top to bottom, the plotted quantities are logarithms of 
density, temperature and pressure divided by the Boltzmann constant ks- In the right 
column, from top to bottom, the panels corresponds to radial velocity, tangential velocity, 
and logarithm of angular momentum. 

Initially, the average gas density in the central part grows rapidly due to the mass depo- 
sition by the stellar cluster. Radiative cooling is enhanced because of its squared dependence 
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on density and thus temperature in the densest zones drops. As soon as temperature de- 
creases to approximately 10^ K, free-bound and bound-bound transitions become the major 
coohng processes. Consequently, when the temperature in some region approaches ~ 3 x 10^ 
K the cooling rate increases steeply and the thermal instability starts to operate. This lowers 
the temperature to 10^ K, and decreases the pressure by three orders of magnitude from that 
(P ~ 3 X 10~^ dyne cm~^) in the hot (~ 10^ K) gas. The cold regions are compressed by 
the surrounding hot gas into dense filaments/clumps. After an initial relaxation period, the 
stagnation radius Rst, which is defined by radiative cooling, remains almost constant. Note 
that in our model this radius determines the amount of gas that flows inwards and that 
later accumulates close to the center and that Rst is different from the Bondi radius which 
is defined by the mass of the central black hole and by the sound speed in the surrounding 
ISM. 

Figure [TJpresents frames of the distribution of the hydrodynamical variables in the whole 
computational domain. The density distribution presents mainly two-phases: hot (few 10^ 
K) gas with low density, and cold (T=10^ K) gas in the densest zones in a form of filaments 
and clumps. The thermal pressure gradient within the stagnation volume {Rgt = 9.2 pc) 
is not high enough to push the cold gas out from the cluster volume and instead the cold 
parcels of gas begin to stream toward the center because of the force of gravity. Due to 
angular momentum conservation, the rotation velocity of the inflowing gas increases as it 
approaches the center, to about 200 km s~^ at ~ 2 pc. Such fast rotation prevents the gas 
to flow further inwards and favors the accumulation of mass around the centrifugal barrier. 
We identify the collection of cold filaments and the dense core at the centrifugal barrier as 
the torus responsible for the obscuration of the central SMBH. While all this happens in the 
central part of the cluster, a stationary cluster wind reaching a terminal speed of about 800 
km s~^ is well sustained above i?st in all simulated cases. Such winds may prevent the inflow 
of matter from larger scales onto the NSC. This suggests that NSC-SMBH may evolve in 
isolation when the feedback from massive stars is highly active. 



4- 1-2. Reference model (model 2) 

The input parameters in this model are the same as in model 1 except that in this case 
we consider the effect of the X-ray radiation from the central SMBH. Figure [2] presents a 
sequence of frames of density, temperature, pressure, radial Ur and tangential U(j, components 
of velocity for the reference model. For comparison, the panels in the right most column 
display the distribution of the corresponding variables at 1 Myr for the case without central 
source of radiation (model 1). As one can note, the stagnation radius remains the same 
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in both models. This imphes that the same amount of mass is accumulated in the central 
part of the cluster. However, the inner zone is affected by the central radiation field, which 
prevents the gas from the fast cooling at the very center (r <1 pc). Therefore in the reference 
model the gas there remains hot (at ~ 10^ K, see the second row in Figure [2]). However, 
the flow at larger radii (1 pc < r < 3 pc) is still dominated by cooling and a substantial 
amount of cold gas is accumulated there forming a filamentary/clumpy torus with a dense 
core supported by rotation at the centrifugal barrier. A fraction of this gas is ionized by the 
SMBH radiation, so, it remains warm (1500 K < T < lO^K). Note that the main difference 
between models 1 and 2 is due to the thermal pressure of the ionized gas in the central zone 
but not due to the radiation pressure. 

Figure |3] shows the multi-phase medium that results from the simulation in the reference 
model. We identify five regions in the log T vs log n plane which represent different compo- 
nents of the flow. Region 1 corresponds to the wind, whose temperature drops due to the 
expansion and radiative cooling. One can also find in this region temperatures that coincide 



with two stable regions (around 10^ and 8 x lO^K) of the Plewa (1995) coohng function 



noticed by Schartmann et al. (2009). Region 2 presents hot (~ lO^K) rarefied gas result- 
ing from shock-shock collisions and radiative heating. Region 3 results from the radiative 
heating of the surfaces of dense clumps in the torus. This region does not exist in the case 
without central source of radiation. Region 4 corresponds to the collection of filaments which 
tend to settle in one of the stable branches of the cooling function. Eventually these cool 
down to the minimum allowed temperature in the simulation due to the increase of density 
as gas moves towards the center. Some high density (n > 10^ cm~^ ) zones are heated up 
to ~ 10^ K by the AGN radiation. These zones are also not present in model 1 which does 
not include the central source of radiation. Region 5 corresponds to the core of the torus at 
the centrifugal barrier where most of the mass is concentrated. The majority of the gas is 
in the cold phase (T<1500 K; see Figure |4]). 



4-1.3. Other models 

In model 3, the stellar cluster is 10 times less massive than that in the reference model. 
Therefore, the density of the inserted gas is an order of magnitude lower. In this case, the 
central source of radiation heats up all the gas within the cluster volume up to few 10'' K 
and prevents the formation of thermally unstable zones (clumps). Therefore the reinserted 
matter does not form torus in this case. 

In model 4, the stellar cluster has the same mass as in the reference model, however 
the cluster is more extended: -Rnsc = 40 pc. Therefore, the density of the inserted mass 
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Fig. 2. — : Time sequence of the distribution of the hydrodynamic variables in the reference 
model. All panels show the z — r plane, with the z axis running horizontally. Each column 
displays frames at the time indicated on the top, for each variable indicated at the end of 
each row. As a comparison, the last column shows the corresponding variables for model 1 
at t = 1 Myr. 

is also lower compared to the reference model. Nevertheless, thermal instabilities occur in 
the densest regions where the accumulated mass forms a torus. However, in this case the 
torus is composed only of warm gas. On the other hand, some clumps formed close to Rst 
eventually join the cluster wind and leave the cluster, reducing the mass accumulation rate. 
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Fig. 3. — : Phase diagram for the standard model. Temperature against number density at 
t = 1 Myr. The diagram was gridded into 70^ cells with the masses of the points depicted 
in color. We identify five components of the flow: 1- wind, 2- hot thermalized gas, 3- heated 
gas in the torus, 4- filaments and clumps, and 5- cold dense core of the torus. The gas tends 
to settle at stable regions in the cooling curve, in particular at around lO^K, however the 
squared dependence of the cooling rate on density leads to lower temperatures. 



In model 5 we consider an extended cluster as in model 4 but, in this case the SMBH is 
10 times more massive and therefore more energetic than in all previous cases. The strong 
radiation field keeps the matter within the whole cluster hot what does not allow to form a 
torus, however a powerful wind as in model 3 is generated. Such cases resemble adiabatic 
calculations given the impact of radiation. 



4.2. Column density and obscuration fraction 

In the AGN models it is supposed that a torus, uniform or clumpy, blocks the light 
coming from the accretion disc. The amount of obscuring gas is usually quantified by the 
column density, i.e. the number of particles per unit area along the line of sight N = J ndl. 
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10^2 cm-2. If > 10^4 cm-2, 



The optical/UV radiation is strongly attenuated above 
the opacity is high enough to block even X-ray photons, in such a case the AGN is said to 



be Compton-thick (Comastri 2004, and references therein). There is observational evidence 
that suggests that a large fraction of AGNs in the local universe are obscured by Compton 



thick gas (Maiolino et al. 1998 Risaliti et al. 1999 Matt 2000) and that most of them are 



associated with Seyfert 2 galaxies (Risaliti et al. 1999). 



Figure |4] displays for the reference model the column density at different evolutionary 
times as a function of the viewing angle as seen from the central SMBH. The line of sight 
along the equator corresponds to ^ = 90°. Different colors correspond to different gaseous 
components. Note that the column density of the warm gas (blue line) is high enough 
(~ lO^^cm"^) to block a large fraction of the X-ray radiation, and thus may turn the AGN 
Compton thick. The cold component (green line) presents gaps and high variability in its 
covering angle, which implies that UV photons can escape from the torus through the holes 
in the neutral gas and an observer can see eventually directly to the center. Such events 
offer a natural explanation to the "mutation" of optically classified Seyfert 2 to Seyfert 1, 



and vice versa. Aretxaga & Terlevich (1994) and Aretxaga et al. (1999) give examples of 



such kind of objects. 

Figure [5] presents the fraction of the sky covered by two different column densities as 
a function of time for models 1, 2 and 4: X24 (red lines) represents the fraction of the sky 
covered by total column densities A^ > 10^^ cm"^, Xcoid,22 (black lines) shows the fraction of 
the sky covered by cold gas with column density A^ > 10^^ cm"^. In the case of model 1, the 
time average values are Xcoid,22 ~ 0.93 and X24 ~ 0.76. In the reference model, Xcoid,22 ~ 0.86 
and X24 ~ 0.92. Note that X24 is larger in the reference model compared to model 1. The 



same tendency was found by Wada (2012) for Eddington ratios in the range 0.01-0.1, where 



more luminous AGNs are obscured over larger solid angles. The fraction Xcoid,22 in the 
reference model is reduced due to ionization by the central source. In the case of model 4 
the torus is composed only of warm gas, with average X24 ~ 0.75 after 1.5 Myr. 



4.3. Comparison with analytic predictions 

The thermally unstable gas inserted within the stagnation radius is attracted by the 
gravity towards the cluster center. Due to its angular momentum it accumulates around the 
centrifugal barrier, where the rotation balance the gravitational attraction. Here we give the 
analytic formula for such radius and show that it is in a good agreement with the numerical 
results. 
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Fig. 4. — : Column densities as seen from the central SMBH for the reference model. The 
red line displays the total column density along each line of sight. The magenta line shows 
the column density for hot gas, T > 3 x 10^ K. The blue line represents the column density 
for warm phase, 1500 K < T < 3 x lO^K. The column densities of cold matter (T < 1500 
K) are shown by green lines. The cold gas does not cover all the sky. 
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Fig. 5. — : Fraction of the sky covered by optically thick gas to optical/UV and X-ray 
radiation as a function of time. Black lines represent the fraction of the sky, Xcoid,22, covered 
by cold gas in the torus with column densities > 10^^ cm~^. Red lines shows the fraction 
of the sky, X24, covered by total column densities > 10^"^ cm~^. 

The radius where the mass accumulates, Rt, is determined by the angular momentum 
of the matter inserted within the stagnation radius Rgt and by the central gravitational 
potential of the SMBH + NSC. Here we neglect the effect of the radiation pressure. The 
value of Rst is defined by radiative cooling, which depends on the mass and compactness of 
the cluster (see TT07 & W08). Thus, in order to estimate the position of the centrifugal 
barrier we need to know Rgt for a given rotation velocity of the star cluster. In the following 
we derive an analytic relation for Rt by assuming a star cluster in solid-body rotation. 

Let us consider the total mass inserted within the stagnation zone at some instant. 
Then the specific angular momentum of a rotating parcel of gas is j = uR"^, where R is the 
projection of r on the equatorial plane, i.e., the distance to the rotation axis. An integration 
over the mass within the stagnation volume V^t = 47ri?g^/3, gives the average specific angular 
momentum inserted within Rgt at some instant: 




(4) 



Then by considering 00 = constant and using spherical coordinates one obtains: 




(5) 



One expects the accumulation of mass at the centrifugal barrier, i.e., at the radius Rt = 
jav/^K, where vk = {GM{Rt) / RtY^'^ corresponds to the Keplerian velocity of the gas or- 
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biting the mass M{Rt) = Msmbh + ^nsc(-Rt)- This leads to the algebraic equation: 

0. 



p4 , ^SMBH p3 p 



The physical solution of equation ^ is: 

.p. 1/2 



a2 t^3 
Jav-'T'NSC 
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(9) 



Thus, Rst is the key parameter to estimate the centrifugal barrier, because it determines 
the amount of mass and angular momentum inserted within the stagnation zone for a given 
set of parameters: -Rnsc? ^nsc? ^smbh, and u. In our calculations, Rst is self-consistently 
determined. It splits the cluster into two zones: the inner one where the reinserted matter 
is accumulated and the outer one where the star cluster wind is formed. It is worth to note 
that Rst is almost independent on the SMBH mass and may have a non-zero value even if 
-^SMBH = 0. This makes our approach different from all the modified Bondi models used so 



Ulrich 


(1976 


); 


Proga & 



Begelman (2003); Krumholz et al. (2005); Inogamov & Sunyaev (2010), where the solution 



depends on the size of the so-called Bondi radius which is a function of the SMBH mass and 
the sound speed in the surrounding ISM. 

Figure [6] shows the analytic predictions for the centrifugal barrier as a function of the 
maximum rotation velocity of the NSC, i.e., the rotation velocity at r = -Rnsc and 6 = 90°. 
Different lines correspond to NSCs of different radii, all of them with a central SMBH of 
10^ Mq. The squares represent the average values in the case of the reference model (-Rnsc = 
10 pc) and model 4 (-Rnsc = 40 pc). In all cases the mass of the NSC corresponds to that in 
the reference model (Mnsc = 3.3 x 10^ Mq). As one can expect Rt increases monotonically 
with the assumed rotation velocity of the cluster. Note that in the range of parameters 
here used, if one considers more extended clusters at a fixed frot, the absolute value of the 
stagnation radius is larger, but the ratio -Rst/-RNSC is smaller. Therefore, the specific angular 
momentum inserted is higher and the mass accumulates at a larger distance from the center. 
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If the NSC is very extended, the impact of coohng is less important and even the absolute 
value of -Rst gets smaller, and eventually the NSC could be in a quasi-adiabatic regim^ In 
such cases Rst tends to a very small value (S04, TT07, W08) and is mainly defined by the 
gravitational potential (S08 and HZIO). 
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Fig. 6. — : Analytic prediction for the centrifugal barrier. The radius Rt where mass 
accumulates as a function of the maximum rotation velocity of star clusters of mass 
Mnsc = 3.3 X IO^Mq and central SMBH of mass Msmbh = IO^Mq, for different star cluster 
radius. The squares represent the reference model (model 2) and model 4, see Table [T| 



4.4. Mass accumulation rate 

In all simulations the hydrodynamic solution reaches a steady state. In the case of our 
reference model it happens after ~0.1 Myr and at about four times longer time in the case 
of 40 pc clusters. From then onwards the matter inserted in the region R^t < r < -Rnsc flows 
through -Rnsc and leaves the cluster as a stationary wind which stops the income of matter 
from a large scale in the galaxy. On the other hand, the mass that remains locked within 




^See, Silich et al. (2004) and Tenorio-Tagle et al. (2007) for a discussion on the threshold energy which 
separates star clusters evolving in the catastrophic cooling regime from those evolving in a quasi-adiabatic 
regime. 
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Table 2:: Mass accumulation rate. 



Model 


^NSC 

(Me) 


^SMBH 

(lO^Mo) 


-RnSC Rst 


^NSC 


A^wind 






(pc) 




(Mo 






1 


3.3 X 10^ 


1 


10 9.2 


31.2 


7.3 


22.1 


1.8 


2 


3.3 X 10*^ 


1 


10 9.2 


31.2 


7.2 


23.9 





4 


3.3 X 10*^ 


1 


40 31.3 


31.2 


15.9 


14.6 






Note. — Mass flow rates for models with typical paranieters of NSC and SMBH in Seyfert type galaxies. 
A substantial amount of the total inserted mass (Mnsc) remains locked within the stagnation radius (i?st), 
and accumulates at a rate ~ Mace in the torus. Note that a good fraction of the total injected mass blows out 
in the NSC wind at rate M^ind- The inflow rate through the inner boundary of the computational domain 
is denoted by Mi,i. 

the stagnation volume streams toward the center and accumulates around the centrifugal 
barrier practically at a constant rate. Table [2] presents the corresponding rates. 

Figure [T] presents for models in Table [2| the absolute and relative (normalized to the 
star cluster mass input rate) values of the rates of mass deposition within the volume of 
the cluster (Mnsc); mass accumulation in the torus (Mace), mass carried away by the wind 
(Mwind); and mass inflow towards the center through the inner zone of the computational 
domain (Mjn). In Model 1: ~ 71% of the total inserted mass accumulates in the torus, 
~ 23% leaves the cluster as a wind. About 6% of the inserted mass escape through the 
inner zone of the computational domain. In the reference model (model 2): ~ 77% of the 
total inserted mass goes into the torus and ~ 23 % goes into the cluster wind. Model 4: 
approximately one half of the inserted mass leaves the cluster as a wind, the rest accumulates 
in the torus. In this case some clumps escape from the cluster, producing peaks in M^ind 
(dashed line) with the corresponding response in Mace (solid line). Note that in models 2 and 
4, radiation pressure prevents the inflow of mass through the inner boundary. The actual 
value of the accretion rate onto the SMBH is beyond the scope of this paper as the inflow of 
gas to the central black hole is inhibited by angular momentum when the viscous processes 



are negligible, Li et al. (2012). Note that in all cases the rate of mass accumulation is given 
by -Rst- 

The mass of the torus at a given time can be estimated from the average Maee- For 
example, at 1 Myr it reaches 2.39 x 10'' M© in the reference model, and 1.46 x 10'' Mq in the 
case of model 4 for the same period. Due to the additional mass loading considered in our 
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models, the mass of the torus grows substantially. Such torus is gravitationally unstable. An 
estimate of the Toomre parameter, Q = fics/SGS, for the obscuring structure in our reference 
model in the region centered at the centrifugal barrier {Rt — 2 pc) with AR = 0.5 pc, sound 
speed Cs —1 km s~^, Keplerian rotation at frequency Q ~ 1.4 x 10~^^ s~^ and surface density 
S = 140 g cm~^ results in Q ~ 5 x 10~^. Therefore, the mass accumulation should lead 
to a continuous star formation, which may amplify the effect of the stellar feedback in the 
nuclear region of the host galaxy. 
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Fig. 7. — : Mass deposition, accumulation and outflow rates. Left and right axis show scales 
of absolute and relative quantities, respectively. Dotted lines represent Mnsc, thick solid 
lines display Mace, dashed lines shows M^ind, and thin solid lines represent Min. The average 
values are given in Table |2j 



5. Impact of the NSC wind on the host galaxy 



In all cases presented in Table [T} the star cluster wind is sufficiently powerful as to 
significantly re-structure the host galaxy ISM leading perhaps to a thick ring along the 
plane of the galaxy, and to a super galactic wind along the host galaxy symmetry axis (as 



m 



Tenorio-Tagle & Munoz-Tunon 1997 1998). 



A simple estimate of the wind power can be obtained from its ram pressure Pram = pw^ 
at the starburst edge. This is in all cases many orders of magnitude larger than the typical 
ISM pressure in our Galaxy (~ 10^^^ dyne cm~^). For example, for models 1 and 2, Pram — 
1.5 X 10"^ dyne cm~^; and about half this value in model 3. In models 4 and 5, Pram = 1-4 
and 2.5x10^^ dyne cm~^, respectively. Such values are comparable with the ram pressure 



Minus /{47tR%sc) 



where ug = [2G(Msmbh + Mnsc)/P: 



NSC 



1/2 



of the freely falling gas, Pg 
is the free fall velocity, in the torus formation model of Hopkins et al. (2012) with the infiow 
rate Mjr, = 10 Mq yr~^ (Pg ~ 2.8 x 10~^ dyne cm~^ for models with 10 pc radius and 
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Pq 8.8 X 10 ® dyne cm ^ for models with -Rnsc = 40 pc). Thus the NSC winds have 
to build up super bubbles and probably super galactic winds preventing, in most cases, the 
falling of the ISM onto the NSC. 



Conclusions 



We present 2D radiation hydrodynamic simulations of the gas reinserted by stars of a 
rotating young NSC with a central SMBH. Our model considers explicitly the impact which 
the stars from a NSC provide on the accretion flow. The model includes constant mass and 
energy deposition from stars and assumes that the mass is inserted with non zero angular 
momentum. It takes into consideration gravity from the central SMBH and from the NSC, 
and accounts for radiative cooling and the heating from a central isotropic source of X-ray 
radiation. 

Here we have shown that the combined effect of gravity from the SMBH+NSC and the 
angular momentum of the inserted mass results in the formation of a dense structure (torus) 
inside the NSC, well within the stagnation radius Rst, defined by radiative cooling. The 
torus is only a few parsecs across, filamentary/clumpy, with a core at the centrifugal barrier. 
It is composed of gas in two phases: a cold phase (T< 1500 K), where dust can survive as 
close as a couple of parsecs from the SMBH, and a warm phase (1500 K<T< 3 x lO^K), 
maintained at this temperature by heating from the central source of radiation. The torus is 
Compton thick and covers a large fraction of the sky, more than 80% in our reference model. 
This obscuring structure is embedded into a low density hot gas. 



Note that models developed by Wada & Norman (2002) and Schartmann et al. (2009) 



and that are here discussed, lead to the flow of the cold gas towards the central zone of the 
host galaxy. The inflow of molecular gas in the inner ~ 150 pc of the Seyfert galaxy NGC 



4051 was detected by Riffel et al. (2008) who suggested that the inflow occurs due to the 



spiral arms that reach the nucleus of the galaxy. Grier et al. (2012) presented the evidence 



for the inflowing gas in the broad line regions of four AGNs: Mrk 335, Mrk 1501, 3C 120 
and PG2130+099. Sim et al. (2012) showed that the parsec-scale inflows do not result in 



significant absorption features in the X-ray spectra since the ionization degree of the infalling 
gas is high. Therefore, the lack of such observations does not rule out our model. 

The accreting mass accumulates in the central region almost at a constant rate, resulting 
after some time in a very massive torus. As soon as it becomes gravitationally unstable a 
second generation of stars may form there leading to the formation of a stellar torus, and 
thus matter would be continuously reprocessed into stars at a rate dictated by the mass 
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accumulation. 

A necessary, but not sufficient, condition for the formation of the torus is that the matter 
reinserted within the NSC evolves in a thermally unstable regime. However, the formation 
of the torus may be prevented by the strong central source of radiation as it is the case in 
our models 3 and 5. 

In all powerful cluster wind is established outside the stagnation radius. Such 

winds can inhibit the income of gas from larger scale in the galaxy. This suggests that during 
the starburst phase, when massive stars dominate the NSC feedback to the host galaxy ISM, 
the NSC-SMBH interplay occurs in isolation. 
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A. The heating rate and acceleration due to the central radiation field 



We consider only the high frequency band of an AGN spectrum, which allows to take 
the advantage of a parametric form of the Compton (Fcompton) and X-ray (Fx) heating 



functions given by Blondin (1994). By means of a ray tracing we calculate the optical depth 
r = JJ^ Kpdr at each radius r, where r-^^ is the inner boundary of the computational domain, 
p is the local gas density, and dr is the radial length of a grid cell. We assume that the 
attenuation in the ionized gas is dominated by Thompson scattering, i.e. the opacity is 
K ~ 0.4 cm^ g~^, and two orders of magnitude higher (see Proga et al. 2000) for the neutral 
gas. 

The optical depth r is used to compute the X-ray flux Fx = Lxe^'^ /{Anr'^), which is 
required to calculate the heating rate and the acceleration due to the radiation pressure. 
The amount of energy emitted in X-rays depends on the SED and the total luminosity of 



the source. For the average AGN SED given in Korista et al. (1997), one gets that about 



8% of the total luminosity is emitted in X-rays. Here we take this fraction and assume that 
the total luminosity corresponds to the Eddington limit for a given SMBH, therefore, we use 
Lx = 0.08LEdd- 
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Once Fx is known, the local acceleration due to the radiation pressure in equation ^ 
is computed as (^rad = Fxi^/c, where c is the speed of light. Then the radial velocity of the 
flow is corrected by graddt. 

The heating rate per unit volume is -^agn = '"■^(rcompton + rx), it is included in equation 
(|3| by increasing the internal energy in each grid cell by HACNdt. Explicitly, rcompton = 
8.9 X 10-3^^(Tx - 4r) and Tx = 1.5 x 10-21^1/47^-1/2(1 _ T/Tx), both in units of erg s"^ 
cm^. Such functions depend on the temperature T of the gas, the characteristic temperature 



Tx of 10 keV X-ray radiation, and on the ionization parameter ^ = AuFx/n (Tarter et al. 



1969). The last parameter characterizes states of ionization equilibrium and depends on the 



local flux and the number density of particles within a grid cell. 
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